Improved Semiclassical Approximation for Bose-Einstein Condensates: Application to 

a BEC in an Optical Potential 
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We present semiclassical descriptions of Bose-Einstein condensates for configurations with spatial 
symmetry, e.g., cylindrical symmetry, and without any symmetry. The description of the cylindrical 
case is quasi-one-dimensional (Q1D), in the sense that one only needs to solve an effective ID 
nonlinear Schrodinger equation, but the solution incorporates correct 3D aspects of the problem. 
The solution in classically allowed regions is matched onto that in classically forbidden regions by a 
connection formula that properly accounts for the nonlinear mean-field interaction. Special cases for 
vortex solutions are treated too. Comparisons of the Q1D solution with full 3D and Thomas-Fermi 
ones are presented. 



PACS numbers: 3.75.FL 03.75.-b, 67.90.+Z, 71.35.Lk 



I. INTRODUCTION 



Simple approximations for describing Bose-Einstein condensates (BECs) have been very useful for understanding 
their physics. For example, in the mean-field approximation, in regions where the local density is large enough, so 
that the mean-field (nonlinear) term in the Gross-Pitaevskii equation (GPE) is much larger than the kinetic-energy 
term, the Thomas-Fermi (TF) approximation offers such a simplified description for the ground state of a BEC in 
a stationary potential [Q. However, in classically forbidden regions of the coordinate space, the density is low, and 
the TF approximation is invalid. It is necessary to match the TF approximation in the region of high density to a 
description valid near the boundaries of the classically allowed motion and in the classically forbidden region for a given 
external potential. For dynamic situations, some simple approximations exist for time-dependent harmonic potentials 
\ P, |j. It would be very useful to have a simple effective one-dimensional (ID) approximation that properly accounts 
for the 3D character of a BEC, for both static and dynamic problems in configurations with spatial symmetry, such 
as a BEC in a cylindrically symmetric potential (e.g., a harmonic trap with cylindrical symmetry, with or without an 
optical potential that varies in space along the symmetry axis of the harmonic potential, see a detailed formulation 
of the model below in Sec. II). Pedri et al. developed a treatment of this kind Q. Another contribution was made 
in Ref. ||, which aimed at a derivation of an effectively one-dimensional (ID) GPE relevant for the description of 
3D BEC by means of the variational approximation. The ID wave function derived in jS) was defined so that it has 
the same density distribution along the symmetry axis of the system as that which could be obtained by integrating 
the distribution produced by the full 3D wave function in the transverse plane. As a result, the effective ID equation 
derived in Ref. B had a non-polynomial nonlinearity. Moreover, the transverse Gaussian distribution adopted in 
Ref. || to derive the ID equation, strongly differs from the actual transverse distribution in the case of large BEC 
' density (which is well approximated by the TF form, see below). 

, Here, our objective is to improve upon the treatments of Refs. |^] in a number of ways. In our treatment in 
k>( ' ^ eC- re gi° ns °f physical space in which the density is sufficient for the application of the mean-field approximation, 
and those in which the density is low, are treated completely differently. In Sec. IV we develop a connection (crossover) 
formula for the wave function in these regions, in analogy with the commonly known formulas between classically 
allowed and forbidden regions in ordinary quantum mechanics ^| (the difference from the connection formula in 
quantum mechanics is due to the fact that the GPE is a nonlinear equation). A WKB approximation for the 
nonlinear Schrodinger equation in an external potential, uniformly valid in classically allowed and forbidden regions, 
was recently proposed in Ref. 0] , but that result does not produce the connection formula. We develop an effectively 
ID treatment of the dynamics in cylindrical symmetry, which properly accounts for the 3D aspects of the problem, 



see Sec. Ill below. In the classically allowed region, the full 3D GPE is reduced to a ID counterpart, the solution 
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to this ID equation being subject to a non-canonical (non-quadratic) normalization condition [see Eq. ( |l3| ) below], 
which is derived from the canonical normalization for the 3D wave function. As a physically relevant application 
of the general method outlined above, in Sec. [v| we calculate the quantum-mechanical distribution of values of the 
longitudinal momentum in a cigar-shaped BEC. This distribution also has a non-canonical form, in comparison with 
ordinary quantum mechanics in ID. In Sec. ^we generalize the approach to treat the case of a BEC with vorticity; in 
this case, the effective ID normalization condition takes on a still more involved form [see Eq. (Q)]. We also develop 
a generalization for the case w hen the cylindrical symmetry is broken by an external potential, so that the effective 



equation is a 2D one, see Sec. VII below. In that case, the 2D wave function is subject to a different non-canonical 



normalization, see Eq. (50) below. In Sec. VIII we provide numerical examples that compare our method with full 
3D calculations for a BEC is static potentials and for dynamics of a BEC in the presence of a time- varying external 
potential. 

We stress that our reduction of the 3D equation to its ID (or 2D, in the broken-symmetry case) counterpart does 
not resort to the Gaussian approximation for the dependence of the wave function on the transverse coordinate(s) in 
the high-density region, where this approximation is not warranted. This is a principal difference from the approach 
developed in Refs. (j, [5j. 

II. A MODEL WITH CYLINDRICAL SYMMETRY 

We begin by defining the kinds of systems of interest to us in the context of a cylindrically symmetric model 
problem. To this end, we consider a BEC in an array of optical traps, in the presence of the gravitational field and 
large-size magnetic trap in the form of a parabolic potential induced by the interaction of the magnetic moment of 
atoms with an external static magnetic field. The static magnetic-trap potential is 

V M (r) = ±m[^z 2 +u;l(x 2 + y 2 )], (1) 

where m is the atomic mass. The optical potential is produced by light beams with identical linear polarizations, whose 
propagation directions lie in one plane with the z axis, forming angles 9/2 and ir — 8/2 with it, (6 — corresponds to 
two beams counter-propagating along the z axis) . Interference between these fields produces a standing- wave potential 
along the z direction, whose amplitude is proportional to the intensity of light. It is assumed that the intensity is 
initially zero and gradually increases with time, hence the light-induced potential experienced by atoms in the BEC 
is 

V L (z,t) = V (t)[l + cos(2kz)}/2, (2) 

where Vo(t) is the optical-potential's amplitude which varies in time, and k = (27r/A p h) sin(#/2) is the wave vector 
of the optical lattices, A p h being the wavelength of light. It is often convenient to discuss the depth of the optical 
potential in units of the recoil energy En — (2nh/\ph) 2 /(2m), which is the kinetic energy gained by an atom when 
it absorbs a photon from the optical lattice. 

Finally, the description of the mean-field dynamics of the condensate is based on the 3D time-dependent GPE, 

ih^(r,t)= \f/2m + V(r,t) + NU \V\ 2 ] *, (3) 
at L J 

where V(r, t) = Vm{y) — mgz + Vl(z, t) is the full potential (the second term is the gravitational potential), p 2 /2m is 
the kinetic-energy operator, Uq — Anaoh 2 /m is the atom-atom interaction strength that is proportional to the s-wave 
scattering length ao, and N is the total number of atoms. Note that, according to Eq. (|l|), the full potential can be 
written as 

V(r , t) = V z (z, t) + V± (fl ) , with V ± (t ± ) = \mu\ (x 2 +y 2 ). (4) 

Equation (|J) can be rewritten, in terms of characteristic diffraction and nonlinear time-scales tjjp and tjvi,) as 
follows [|[ g: 



~dt 



r ^V 2 -V(r,t)/K-^P 

foF *NL 



* . (5) 



Here, the diffraction time toF — 2mr\ F /h, with ttf = yM^I '(rnui 2 ) and cj — (lJzLU 2 ^) 1 / 3 , and the nonlinear time 
£nl = [G\^ m \ 2 /Ti)^ 1 = (fi/h)^ 1 can be expressed in terms of the chemical potential, where \^ m \ is the maximum 
magnitude of W Another useful length scale is ttf.z = \/2[i/ (muSf), i.e., the TF radius in the z direction. This is 
the size of a TF wave function in the z direction for the harmonic potential muj 2 z 2 /2. 
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III. THE SEMICLASSICAL APPROXIMATION 



We consider a semiclassical approach based on the Thomas-Fermi approximation for a 3D BEC wave function in a 
cylindrically symmetric potential. Our treatment is broken up into different approaches depending upon whether the 
atomic gas density is high (in the classically allowed regions of the coordinate space) or low (in classically forbidden 
regions). 



A. Classically Allowed Region 

In the classically allowed region [not too close to its boundaries so that the atomic gas density, and hence, the 
nonlinear term in Eq. (||), remain sufficiently large], our description is based on the following ansatz for the 3D wave 
function *&(]*, t), 

n,, t )=*^)( G - v ^f^ t)? f (o) 

where ip(z, t) is a newly defined effective ID wave function, and G = NUq [cf. Eq. (g)]. Note that the ansatz assumes 
a fairly simple relation between the squared ID and 3D wave functions, 

M 2 = Vl/G + |*| 3 . (7) 

Of course, the ansatz ^ makes sense in the region where 

G\^(z,t)\ 2 -V ± (r ± ) >0. (8) 

Substituting the ansatz (Q) into the 3D GPE and neglecting the transverse part of the kinetic-energy operator in 
the spirit of the TF approximation, we arrive at an effective ID GPE, 



dip 



^-fi+V x (z,t) + G\rl>\ 2 
2m 



4>, (9) 



upon neglecting terms proportional to V±. Strictly speaking, to derive Eq. ([)]) we need a condition V± -C G^ 2 , in 
terms of the 3D wave function. However, it will be shown below that the ID equation (^) is a reasonable approximate 
model even when the terms V± and C| \I/ 1 2 are on the same order of magnitude. The full TF approximation for 
stationary solutions to Eq. (^) can be obtained, as discussed in Sec. Ill C below, if one further neglects the longitudinal 
kinetic-energy operator (l/2m)p 2 ,. We stress that ansatz (0) is not an exact solution of the full 3D GPE, and it cannot 
describe radial excitations of the BEC. 

For the case of a harmonic transverse potential V± from Eq. (Q), the above condition (||) yields 

2C 

r 2 <r 2 m (z,t) = \^( z ,t)\ 2 , (10) 
muj ^ 

where r is the radial coordinate in the plane (x,y). The quantity r m introduced in Eq. ([To| ) may be regarded as a 
definition of the radius of the cigar-shaped BEC at a given values of z and t. The 3D wave function \I>(r, t) is subject 
to the ordinary normalization condition, 

2tt / rdr I dz |*(r, z, t)\ 2 = 1. (11) 

J0 J-oo 



Upon noting that the integration over r in Eq. (|llj) is confined to the region r < r rn and substitution of the ansatz 
(||) into Eq. (|ll|), we obtain a result, 

/+oo rr m r+°° r r ™ 

dz rdr\f(r,z,t)\ 2 = — dz rdr [G\ip(z, t)\ 2 - muj'ir 2 /2] 



dz l rdr (r 2 m - r 2 ) = ( + °° \^(z, t)\ A dz = 1. (12) 



G / M L > mw 2 
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Thus, according to Eq. (J12|), the usual normalization condition for the 3D wave function, Eq. (|llj), generates the 
following non-canonical normalization condition for the effective ID wave function: 

r +QO mu 2 

/ \^t)\Uz = ^. (13) 

We stress that this abnormal-looking condition is a direct result of the standard full normalization condition ( p"l| ) and 
the ansatz (^) adopted for the 3D wave function. 

We note that, as follows from Eq. (||), ip(z,t) = = 0, y = 0, i.e., the function ip is the particular value of 
the full wave function 'J on the axis x = y = (therefore, the functions ip and VP are measured in the same units). 
Despite the on-axis identity between the functions \I/ and ip, the latter one does not have the interpretation as a 
probability amplitude for the distribution of atoms along the z axis; instead, the probability for finding a particle in 
the region between z and z + dz (integrated in the transverse plane) is 

POO s~1 

P(z)dz = 2ir rdr\^(r,z,t)\ 2 dz = T \ip(z, t)\ 4 dz , (14) 

Jo muj ± 

cf. Eqs. ([|) and (|l|). 

Recall that, in contrast to linear quantum mechanics, the normalization of the wave function ip is important and 
affects physical results in nonlinear theories of the GPE type. Indeed, the strength of the nonlinear mean-field term 
in Eq. (^|) is determined by the maximum value of \ip\ 2 and is thus affected by the normalization. 



B. Stationary and Slowly Varying Cases 

A stationary solution to the 3D GPE with a time independent potential V(r) can be approximated, in the classically 
allowed regions, by a stationary version of the ansatz (g), ip(r, z, t) = <p(z) exp [— (i/h) fit], where the function <p{z) 
satisfies the equation 



(15) 



following from Eq. (||). 

For problems with a slow time variation, we can consider an instantaneous eigenstate of the nonlinear time- 
dependent GPE equation. Adiabatically varying potentials V z {z,t) can be treated by calculating the instantaneous 
chemical potential and quasi-stationary wave function (p{z; t) in the instantaneous external potential, and then forming 
the full time-dependent solution as 



-0(r, z, t) — <p(z] t) exp 



(i/h) / dt'ix{t') 



(16) 



Strictly speaking, the non-canonical normalization condition (|13|), unlike the canonical (quadratic) one, is not com- 
patible with the full time-dependent effective ID GPE (||). However, there is no problem with the compatibility in 
the case of the adiabatically slow evolution. 



C. Full Thomas- Fermi Approximation in the Classically Allowed Region 



The full 3D TF approximation can be recovered if we apply the TF approximation directly to the ID equation 
neglecting the kinetic-energy operator in it, so that the solution will be 



iP(z,t) 



H-V g (z,t) 



G 



exp 



i 

h Jo 



dt'n(t') 



According to Eq. (^|), this yields the full 3D TF wave function, 

^-V ± (r) ± -V z (z,t)\ 1/2 



*M) 



G 



exp 



dt'n(t') 



(17) 



(18) 
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Substituting Eq. ( |l7| ) into the non-canonical normalization condition (|T^) yields 

[fi(t) ~ V z {z,i)fdz = Gmul/n, (19) 



where the region of integration over z is restricted by the condition fi(t) — V z (z, t) > 0. Hence, the normalization 
condition (|l^) determines the chemical potential n(t). Note that the condition ( |l9| ) is equivalent to the usual form of 
the condition which determines the chemical potential in the framework of the TF approximation applied to the full 
3D equation (ph: 



rdr / dz \p(t) - V(r,t)] = 1, (20) 



2tt ' x 
~G 

where the integration is performed over the region in which ji{t) — V(r,t) > 0. Performing the algebra, we arrive at 
the usual expression for the chemical potential in the static harmonic 3D potential (without an optical component): 

M=i[15G/(47r)] 2/5 (m^ 2 ) 3 / 5 . 

Finally, for any potential V z (z) (and the harmonic potential V±), the effective probability density defined by Eq. jl^ ) 
takes the following form in the TF approximation: 



p tf(*) = — ^|m-^(z)I 2 , (21) 

in the region where fi — V z (z) > 0; otherwise, Ptf(z) = 0. 



D. Classically Forbidden Regions 

In the classically forbidden regions, the density of atoms is small, therefore the nonlinear term in the GPE may be 
dropped, so that it becomes tantamount to the ordinary quantum-mechanical Schrodinger equation, hence we adopt 
the following product ansatz for $>(r,t): 

*(r, t) = iP(z, t) exp(-ri/2i?l - iw x t/2) , (22) 

where the transverse squared radius is R\ = fr/muj±. We stress that the Gaussian approximation for the transverse 
part of the ansatz ( p2| ) is appropriate, unlike in the classically allowed region, as the equation is effectively linear in 
the present case. Upon substituting Eq. ([22]) into the linearized GPE, it is straightforward to obtain an effective ID 
linear Schrodinger equation, 



l-p 2 z + V z (z,t) 



4> ■ (23) 



Quite naturally, Eq. ( J23| ) is equivalent to Eq. (M) in the classically forbidden region, as in this region the nonlinear 
term in Eq. (||) is negligible. 

IV. MATCHING CONDITIONS 

Equations (^|) or (^) can be rewritten in the following form by rescaling t, z and tp: 

1 2 
iu t + -u zz - U(z, t)u - 7 \u\ u = 0, (24) 

where 7 > is a properly normalized nonlinear coupling strength. A solution to Eq. ( p4| ) is sought for as 

u(z,t) = v(z;t)enp J dt'fx(t')^ , (25) 
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with a real chemical potential /j,(t) and a real function v(z;t), cf. Eq. (|16|). Here, the time dependence of U{z,t) is 
presumed to be slow enough, and v(z; t) satisfies a quasi-stationary equation, 



Mt)-U(z,t)]v + ^-"fv 3 



= 0. 



(26) 



As in the ordinary semiclassical form of quantum mechanics, it is necessary to match the approximations for the 
wave function across the classical turning point, which separates the classically allowed and forbidden regions in the 
ID space. As well as in linear quantum mechanics, the wave function in the latter region will be taken in the WKB 
approximation, see Eq. (^7|) below. However, a crucial difference from the standard theory is that the wave function 
in the classically allowed area is taken not in the corresponding version of the WKB approximation, but rather in the 
TF form. This, of course, drastically changes the matching problem (see also Ref. @). 

Deeply under the barrier, i.e., for large positive values of the potential U(z), the density of particles is small, hence, 
as it was already mentioned above, the nonlinear term in Eq. (|26| ) may be dropped, and a solution may be presented 
in the standard semi-classical (WKB) approximation, 



v(z) 



C 



[2(C/(z)- M )] 1/4 



exp 



a/2 (U(z') - n) dz' 



(27) 



where, for the defmiteness, we choose zq as the classical turning point, at which U(zq) = fi, and C is (for the time 
being) an arbitrary real constant (because C is arbitrary, zo may indeed be chosen arbitrarily). It is also assumed 
that the classically forbidden region is located at z > Zq, i.e., to the right of the turning point. As usual, the WKB 
approximation (Ej) is not valid too close to the turning point. 

On the other hand, the solution in the classically allowed region (z < zo), not too close to the turning point, is 
taken in the usual TF approximation as described above, cf. Eq. (17): 



'V-U(z) 
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(28) 



Note that, unlike the WKB solution (p7f), the TF approximate solution (28|) does not contain any arbitrary constant. 

Now we need to match the two approximations (27) and ( p8| ) across the turning point, in a vicinity of which both 
approximations are not applicable, the eventual objective being to find the constant C in Eq. (|27|). Following the usual 
quantum-mechanical approach, one can cast the matching problem into a standard form, expanding the potential in 
a vicinity of the turning point, 



/j - U(z) w F (z - z ) 



(29) 



where Fq is the value of the potential force at the turning point (in the present case, Fq < 0). The accordingly 
modified version of Eq. (^6|) is 



dz 



which is transformed into a normalized form, 



—j + F a (z - z ) v - 2 7 t> 3 = 0, 



by means of rescalings 



d 2 w 

~d~e 



^(2|F |) 1/3 (z 



£w + 2w 3 , 



zo 



_ (2]Jbj) 1/3 ; 
V7 



(30) 



(31) 



Equation ( pl| ) is a particular case of a classical equation known as the Painleve transcendental of the second type. 
The full form of this equation (a standard notation for which is Pn) is 



d 2 w 



-j^- = £,w + 2w 3 + a, 



cf. Eq. fl30P, where a is an arbitrary real parameter; in the present case, a = 0. The use of the expansion (|29j) 
and simplified equation (BO) for matching different asymptotic solutions of an equation equivalent to Eq. (EQ) was 
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proposed, in a context different from BEC, in Ref. jTo) (however, the asymptotic form of the solution in the classically 
allowed region, for which the analysis was done in Rcf. [|o| , was different from Eq. (|2|) : it corresponded to a nonlinear 
wave function oscillating in space, rather than to the TF case). Here, it is necessary to find a solution to Eq. ( |30| ) 
with the property that it takes the asymptotic forms 

w « and w « ^ exp (- ^ 3/2 ) (32) 

at £ — > — oo and £ — > +00, respectively, in accordance with Eqs. (p8h and (|27|). The actual problem is to find a value 



of the universal numerical constant C in Eqs. (32) from matching the solution to its uniquely defined asymptotic form 
at £ — * —00. 

As mentioned above, this matching problem is different from its counterpart in ordinary (linear) quantum mechanics. 
An exact solution to this problem is available in the mathematical literature (see Ref. |0| and references to original 
works therein) . The final result of the analysis is 

C = — ^ sa 0. 282 . (33) 
2y^ 

Using this exact result, and undoing the rescalings (|3l|), we obtain the value of the constant C in the WKB solution 



C = 5^Ua/& (34) 
With the relation p4), the WKB expression (E7]) for the wave function under the barrier is completely defined. Of 




course, the result presented in the form (34) also applies when the classically forbidden region is located to the left 
(rather than to the right) of the classical turning point. 

This result for the matching problem applies as well to the case of the adiabatically slow variation of the parameters, 
e.g., when Fq slowly varies as a function of time (or when 7 changes with time due to variation of the trap potential 
with time). The result does not apply to the case of a vertical potential wall (when, formally, Fq = 00). However, 
in this case, the solution is almost trivial. Indeed, assume that at the point z — there is a jump of the potential 
from a large negative value U- inside the well (z < 0) to a large positive value U+ in the classically forbidden region 
(z > 0). Then, in the allowed region, the TF solution in the form v = y/JJi — U- ) / 7 is valid everywhere up to the 
turning point, the exact solution in the forbidden region is 



v = C exp {-y/2{U+ - Jijz) 



and the constant C is immediately found from the continuity condition, C = y/ (ji — U-)/j. 



V. MOMENTUM DISTRIBUTION 



The approach developed above can be naturally applied to calculate the distribution of values of the longitudinal 
momentum p in a given quantum state |^/), which can be measured in a direct experiment. The distribution is 
determined by the scalar product 

V(p) = \(m\\ (35) 

where, up to a normalization factor, 

(p|=exp(=jp) (36) 

is the conjugate eigenfunction of the momentum operator. Thus, according to Eq. (|35|), to determine the momentum 
distribution in the z direction we need to calculate nothing else but the ID Fourier transform of the 3D wave function 
(0), additionally integrated in the transverse direction [the scalar product in Eq. ( |35| ) assumes, of course, the full 3D 
integration] : 

(p\f)=2irj dzexp(^j^J J *(rj_, z,t) r ± dr ± , (37) 
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where the multiplier 2ir is generated by the angular integration in the transverse plane. 

Substituting the expression @ for \P in the classically allowed region into ([37|), the integration over r± is confined 
to the interval < r± < r m , where r^(z,i) is the same as defined above by Eq. (10). Taking into account the form 
of the transverse potential (^) , we arrive at an expression 

= 2n J °° r 2 m dz exp (^j^ *Kz, t) j ^/T^pdp , (38) 

where p = r^/r m . Further, using an elementary formula — p 2 pdp = 1/3 and the expression (10) for r m , we 

obtain from Eq. (p8|), 

(p|*) = ^-^r f +OC dzexp(^-)i>(z,t)\i>(z,t)\ 2 , (39) 



3 mJ\_ J _ 00 



which is the final result: E q. (|39| ) tells one that the amplitude (pl^) of the probability distribution, that determines 
the probability as per Eq. (pq), is 



(pI*) = -^ — T-^iv-rxp/n), (40) 



where J 7 is the symbol of the Fourier transform, 



1 



+ 00 



Hf}<y>) = j^yji J ex P H<^) • ( 41 ) 

In the classically forbidden region, the usual quantum mechanical expressions for the momentum distribution apply. 

Expectation values of operators involving only the longitudinal momentum can be computed as follows, if we neglect 
a contribution from the classically forbidden regions. Suppose we have an operator O(p), such as, for instance, p itself 
or the kinetic energy, p 2 / (2m). Then, the expectation value of the operator 0{p) is given by 

/+oo 
dpO(p)\(p\^)\ 2 . (42) 
-oo 

The expression (|40|) should be then substituted for (p\^>) in Eq. (Eg). 



VI. THE CONDENSATE WITH VORTICITY 



The above consideration can be generalized for a case when the condensate inside the cylindrically symmetric region 
is given vorticity, so that the wave function has the form 

*(r, z, t) = exp (iW) $(r, z, t), (43) 

where 9 is the angular coordinate in the transverse plane, and I is the vorticity quantum number. The ansatz @ is 
then replaced by 



¥(r, z, t) = exp (iW) G-^ 2 ^G\ilJi(z,t)\ 2 - (muir 2 /2 + 



h 2 P 

2mr 2 



The corresponding normalization condition replacing Eq. (O) takes the form 



2n 



dz 



rdr |^(r,z,t)| 2 = 1, 



where now, 



raw | 



G\M 2 + ^G 2 \M 4 -{l^iS' 



G\M 2 - V G2 I^! 4 -(lhu±Y 



(44) 

(45) 

(46) 
(47) 
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The corresponding ID GPE equation for ipi(z,t) is 



ih 



dt 



2m 



p 2 z + V z (z,t) + G\^ 



Finally, the substitution of the expressions (|44j ) and fl46|), ( |47| ) into Eq. ( [45| ) yields, after a straightforward calculation 
of the integral over the radial variable r, an effectively ID normalization condition in a complicated form, which is a 
generalization of the above non-canonical normalization condition ( |l3| ) corresponding to I = 0: 



-foo 



dz 

oo 



G y/ G 2 |^)! 4 -(^) 2 

2 4 



\muj_ 



In 



2G 2 |^)| 4 l | 2G\Mztf G*\Mz)\ A 



- 1 



2G 



(48) 



Note that, with I — 0, this result indeed reduces to Eq. (13). A consequence of Eq. 



take values smaller than 



\ / 1 



is that cannot 



(lhiv±) /G. Furthermore, from Eq. it is clear that the condition 



G\ipi(z,t)\ 2 - muj 2 ± r 2 /2 + (U/r) 2 / (2m) 



> must be satisfied for the ansatz 



to be valid. 



VII. THE BROKEN-SYMMETRY CONDENSATE 



As mentioned in the introduction, the approach developed above can be extended to the case when the cylindrical 
symmetry about the z axis is broken by a non-axisymmetric potential V xy (x,y). In this case, essentially the same 
ansatz as given by Eq. may be employed if the potential is harmonic in z. Accordingly, we take the transverse 
potential as V±(z) = (1/2) rauP^z 2 . cf. Eq. (^), and an ansatz in the form 



*(r,t) = ip(x,y,t) 



G -muj 2 _z 2 /2\%l)(x,y,t)\' 
G 



1/2 



(49) 



Note that the interpretation of the function ip(x,y,t) is similar to that of the function -0(z,i) in the ansatz (Q): it 
coincides with the full 3D wave function ^(x, y, z,t) at z = 0. Substitution of Eq. J49| ) into the 3D normalization 
condition (nW and straightforward integration in the z-direction yields a result 



\i/j(x,y)\ dxdy = 



2G 



(50) 



cf. Eq. (|l^) . This is another example of the non-canonical normalization condition which is generated by the reduction 
in the effective space dimension. As for the effective 2D GPE generated by the ansatz (0), it has the usual form, 



dip(x,y,t) 
dt 



1 

2m 



(pI+pI) + v x 



(51) 



cf. Eq. (||) [a formal condition under which Eq. ( ^l| ) can be derived from the 3D equation (Q) is muP^z 2 <§; \ip(x, y, t)\ 2 ]. 

We do not consider the problem of matching the wave function in the classically allowed and forbidden regions 
in the framework of the 2D equation (pl|), as the WKB approximation for the classically forbidden region is itself 
problematic in the 2D case. In fact, this approximation was only elaborated for the 2D motion in an axially symmetric 
field which does not correspond to the situation of interest in the present context [the axial symmetry in the (x, y) 
plane will be destroyed by the optical lattice] . 



VIII. NUMERICAL RESULTS 



For the numerical solution of Eqs 
putational grids had 65536 and 32 x 



(|^) and (^|), we used the standard split-step operator method | 
32 x 2048 points for the ID and 3D cases, respectively wit 



|. The corn- 
spatial steps 



Ah z ~ 0.0002 7"tf z for the ID geometry and Ah x — Ah y w 0.13 ttf and Ah z « 0.005 ttf for the 3D case. We used 
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two different methods for finding stationary solutions to GPE ([la), and to its 3D counterpart: the hrst technique 
was an imaginary-time version of the split-step operator method [Q , and the second is the standard finite-difference 
method used to solve a two-point boundary value problem [jl3|. It is important to note that, by treating Eq. ( |l5| ) 
as a two-point boundary problem, the initial conditions must include the eigenvalue fi. This will, generally, give an 
unphysical solution for <fi since the normalization condition will not be met. By gradually changing the value of /x, and 
redetermining the wave function of the nonlinear GPE using the finite-difference method, we can follow the surface 
of solutions to Eq. ( |l5| ) until a physical solution is obtained. 

Numerically, we find the imaginary time split-step relaxation method painfully difficult to converge in 3D for the 
cigar-shaped geometry that we used in the calculations, both with and without the optical potential. The finite- 
difference two-point boundary value method appears to be more efficient. For dynamical simulations, however, we 
used only the split-step method. 

We consider a 87 Rb condensate with 10 6 atoms in a harmonic potential with trap frequencies lu z = 100 Hz and 
lj x = Ljy = 20 Hz, and s-wave scattering length ao = 5.017 nm. We first consider the case when only the harmonic 
potential is present; the action of an optical potential on the trapped BEC will be considered below. For the 
parameters used in the calculations, the TF radius in the z direction is ttf,z — 36.35 /an, and the chemical potential 
is fj, = 1.507 x 10~ 30 J. 

Figure |l| shows the probability P(z) versus z, as calculated using Eq. (@). Also shown is the probability Pyf(z) 
as obtained from the full TF approximation based on Eq. ( |2~L| ), and a probability Pzd(z) found from the numerical 
solution of the full 3D GPE. The curve obtained from the 3D GPE lies on top of the TF curve, being nearly 
undiscernible from it (it is no surprise that TF is a good approximation for 10 6 Rb atoms in a harmonic trap). 
Clearly, the comparison of the P(z) found by means of the approximation developed above with the TF and full 3D 
results is excellent. A slight deterioration occurs in the region where the density is low. 

Figure ^ shows the probability distribution for the same harmonic potential as in Fig. [l], but with an added optical- 
lattice potential as given by Eq. (^), with the wavelength A p h = 840 nm, relative angle 6 = 10 degrees between the two 
light beams, and the constant amplitude Vo = 5Er. For this wavelength, the recoil energy Er = ft 2 (27r / A) 2 / (2m) = 
2.15 x 10~ 30 J, so the optical-potential's strength is about 1.4 times the value of the chemical potential in the absence 
of the optical lattice. In the presence of the optical lattice, the chemical potential is calculated to be 3.15 x 10 -30 
J (slightly more than twice the chemical potential without the optical lattice) . Also shown in Fig. || are the optical 
potential Vl(z), and P(z) determined without the optical lattice (as in Fig. |l|). As is seen from the figure, the optical 
potential squeezes the atomic density out to larger z\ it also squeezes it out to larger x and y, see below. 

Figure || shows the comparison of P(z) with Ptf(z) and P3d(z), where Ptf(z) is calculated using Eq. (21), and 
with P3d(z) as obtained by solving the stationary 3D GPE. The TF result still provides an adequate description, 
despite the fact that the kinetic energy is more important in this case than in the case without the optical potential. 
Hence we conclude that the kinetic energy remains much less significant than the potential and the mean-field energy 
in this case. Our method produces results closer to the 3D result than the TF. The kinetic energy does broaden 
the wave function in each optical well, hence peaks of the wave functions are lower than in the TF approximation, 
as evident in Fig. ^. Note, however, that our method, being based on a TF-type approximation in the transverse 
dimension, restricts the diffusion of the wave function due to the kinetic energy to the z direction (but the wave 
function is definitely squeezed into the transverse dimension due to the optical potential and the mean-field - see next 
paragraph). The inset in the figure is a blowup of the region near z/rxF.z = 1- 

Figure f| shows |^(r)| 2 versus r and z. A striking aspect of this figure is the extent to which the wave function is 
squeezed out to larger r in the presence of the optical potential. Without this potential, the size of the wave function 
in r is ttf.x (= ^tf,j/), but now it is squeezed out to about 5rxF,:r- The size of the wave function in the radial 
direction depends upon z, as does the extent of the squeezing in r. The distribution of ^(r)! 2 versus r and z, as 
produced by the stationary 3D GPE, is similar to that obtained by our method. 

Lastly, in Fig. [| we show results of a dynamical calculation in which we varied the optical potential as a function 
of time, so that the peak strength depended on time as Vo(i) = 5Er x cxp[— ((t — tp)/a) 2 ], with tp = 1.12 ms and 
a = tp/2 = 0.56 ms. The calculated probability distribution P(z,t) is shown at four different times, viz., at t = 
before the optical potential is ramped up, at t = tp/5 when the optical potential is still rather small, at t = tp/2 
when the optical potential is somewhat less than 0.4 times its final value of 5Er, and at t = tp, when Vo(tp) = 5Er. 
The dynamics are initially adiabatic, but, clearly, by the final time, tp = 1.12 ms, the dynamics cease to be adiabatic 
(compare the result with the probability distribution shown in Fig. 0). The dynamics of the probability P(z, t) versus 
z, calculated using the 3D GPE, are similar to that shown in Fig. u. 
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IX. SUMMARY AND CONCLUSIONS 

We have proposed several improvements to the semiclassical description of BECs in three dimensions. First, an 
ansatz that makes it possible to reduce the corresponding 3D GPE to an effectively ID equation in the cylindrically 
symmetric case was put forward. An interesting feature of this approach is that the corresponding ID normalization 
condition, which follows from the standard normalization condition in the full 3D description, takes a non-canonical 
form, containing the fourth power of the ID wave function, rather than its square. Also non-canonical is an expression 
for the probability density of the distribution of the ID momentum. These results were further extended to cases when 
the BEC has vorticity, and when cylindrical symmetry is absent; these cases yield additional examples of non-canonical 
normalization conditions, sometimes of quite complicated form. 

Another result, obtained in the framework of the effectively ID description, is an explicit matching formula between 
the TF approximation valid in the classically allowed region, and the exponentially vanishing WKB approximation 
valid in the classically forbidden region. Here, an exact solution to a problem, found long ago in an abstract mathe- 
matical context, determines the arbitrary constant in front of the exponentially decaying WKB wave function. 

To verify the validity of analytical approximations developed in this work, we have performed direct numerical cal- 
culations of bound states and of dynamics in a time-depenendent potential, and compared the probability distributions 
obtained with full 3D results. The comparison shows that the analytical approximations are quite accurate. 
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FIG. 1: The probability P(z) for the 3D harmonic potential versus z/r tf,z, as found numerically using Eq. (g), and Ptf(z) 
versus z/rTP, z , as calculated using Eq. (^). The probability P'jd(z) found from the stationary version of the full 3D equation 
(see the text) cannot be discerned, as it lies on top of Ptf(z). 
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FIG. 2: The probability P(z) for the combination of the 3D harmonic and optical potentials vs z/vtf,z, as found numerically 
using Eq. (^). For comparison, the probability P{z) without the optical potential, and the optical potential itself are also 
shown. 
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FIG. 3: The probability distributions P(z), P3d{z) and Ptf(z) versus z/Ytf,z for the combination of the 3D harmonic and 
optical potentials. The inset is a blowup of the region near z/Vtf.z = 1. 
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FIG. 4: The probability |*(r)| 2 vs x and z. 
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FIG. 5: The probability P(z,t) vs z at four different values of time, t = [Vo(0) = 0], t = t F /5 [V (t F /5) = 0.387 E R ], t = t F /2 
[Vo(t F /2) = 1.839 t = t F [Vo(t F ) = 5E R ]. 



